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1.0  INTRODUCTION 


1.1  Aims  of  the  Research 

Solid  hydrogen  doped  with  atomic  impuritiesh2  (e.g.,  Li,  B,  H)  is  potential  high  energy 
density  matter  (HEDM)  for  use  in  rocket  propulsion.  One  critical  issue  in  these  materials  is  the 
metastable  structure  of  the  impurity  trapping  sites.  This  issue  was  a  primary  focus  of  research  in 
the  HEDM  Program  over  the  past  five  years  (see  Sec.  2.0  Research  Results,  Data  Obtained  and 
Accomplishments).  However,  a  second  and  probably  even  more  important  issue  is  the  rate  of 
diffusion  and  recombination  of  impurities  in  the  solid  hydrogen  host.  This  second  issue  became  the 
focus  of  the  research  in  the  final  two  years,  and  it  is  presently  the  focus  of  the  Air  Force  Office  of 
Scientific  Research  (AFOSR)-funded  continuation  of  the  project.  Indeed,  the  ability  to  "scale-up" 
potential  cryogenic  HEDMs  in  order  to  achieve  an  impurity  concentration  of  a  few  mole  percent  - 
and  to  have  this  material  be  stable  for  a  reasonable  period  of  time  -  is  perhaps  the  central  issue  that 
must  be  addressed  before  such  materials  can  be  usefully  employed  in  propulsion  systems.  Thus, 
the  timescale  for  recombination  of  the  energy-storing  impurity  species  within  the  condensed  phase 
host  is  the  key  physical  parameter  that  characterizes  the  overall  stability  of  the  HEDM.  Ideally,  the 
timescale  for  the  recombination  reactions  should  be  at  least  on  the  order  of  days,  depending  of 
course  on  the  specific  use  of  the  material. 

Theoretical  studies  can  help  us  to  understand  the  key  features  that  govern  the  stability  of  a 
HEDM.  For  example,  the  average  rate  of  impurity  recombination  in  condensed  matter  is  usually 
governed  by  two  dynamical  timescales:  the  rate  of  impurity  self-diffusion  and  the  intrinsic  (i.e., 
local)  rate  of  impurity  recombination.  The  overall  recombination  process  in  low  temperature  solid 
HEDM  may,  therefore,  be  characterized  as  a  diffusion-influenced  chemical  reaction  having  a  rate 
constant  for  the  macroscopic  recombination  rate.  This  fundamentally  important  constant  -  which 
has  become  the  target  of  the  theoretical  investigations  -  is  approximately  given  by  the  equation^-^ 

(1) 

where  is  the  intrinsic  recombination  rate  for  the  two  impurities  (i.e.,  after  reaching  some 
"contact"  distance  rj  and  is  the  rate  at  which  the  two  impurities  diffuse  into  proximity  with  one 
another.  The  latter  rate  is  related  to  the  impurity  self-diffusion  constant  D  such  that  k^  ~  Akt^D. 

As  may  be  readily  verified  by  inspection  of  Equation  1 ,  the  impurity  self-diffusion  process  limits 
the  rate  when  kj)  «  ki„,  whereas  the  microscopic  recombination  step  is  rate  limiting  if  «  kjj. 

Both  the  self-diffusion  and  recombination  dynamics  of  atomic  impurities  in  solid  matrices 
involve  complicated  dynamical  processes.  In  such  processes,  a  particle  may  need  to  surmount  (or 
tunnel  through)  a  barrier  in  the  effective  potential  energy  surface  along  some  "reaction  coordinate" 
in  order  to  reach  a  new  stable/metastable  state.2-i3  For  the  diffusion  process,  the  guest  atom  may 
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move  in  a  series  of  uncorrelated  hops  between  the  potential  energy  minima  in  the  condensed  phase 
host,  or  it  may  follow  a  mobile  vacancy.  In  the  recombination  step,  the  diffusing  atom  may 
undergo  an  activated  transition  to  occupy  the  same  site  as  a  second  guest  atom,  allowing  a  bond 
to  be  formed  with  that  atom.  The  latter  process  would  be  quite  exothermic  in  a  cryogenic  HEDM. 
The  interaction  of  the  diffusing  particle  with  the  condensed  phase  matrix  provides  an  overall 
macroscopic  free  energy  barrier  to  recombination,  thereby  leading  to  a  metastable  material. 

Based  on  the  above  picture,  it  became  apparent  that  a  better  microscopic  understanding  of 
both  the  impurity  diffusion  and  recombination  rates  in  low  temperature  HEDM  was  (and  still  is) 
necessary  in  order  to  characterize  their  stability.  The  present  contract  was  the  first  condensed  matter 
theoretical  project  funded  by  the  Air  Force  to  study  the  behavior  of  cryogenic  HEDMs  when  this 
problem  became  a  research  priority  in  the  early  1990s.  Expertise  in  condensed  matter  theory  and 
simulation  was  brought  to  bear  on  the  problem,  as  well  as  expertise  in  high  performance 
computing,  which  proved  to  be  an  asset  to  the  HEDM  Program.  During  the  period  of  the  contract, 
the  focus  was  to  develop  and  apply  both  quantum  path  integral  simulation  and  analytical  methods, 
with  the  early  goal  being  the  modeling  of  the  trapping  sites  and  spectroscopy  of  reactive  impurity 
species  such  as  lithium  or  hydrogen  atoms  in  the  solid  para-hydrogen  matrix.  These  efforts  were  in 
no  small  part  motivated  by  the  experimental  work  of  Fajardo ^>2  and  co-workers  at  Phillips 
Laboratory  at  Edwards  AFB.  The  computational  studies  in  the  second  part  of  the  project  were  then 
aimed  at  understanding  the  microscopic  factors  that  directly  influence  the  fundamental  dynamical 
process  of  impurity  diffusion  and  recombination.  These  studies  required  a  novel  theoretical/ 
computational  methodology  to  carry  them  out,  as  explained  in  the  next  paragraph. 

The  cryogenic  HEDM  systems  based  on  solid  hydrogen  (e.g.,  Li  in  H2)  involve  atomic 
impurities  embedded  in  a  low  temperature  solid  hydrogen  matrix.^^’^^  Under  normal  conditions 
hydrogen  solidifies  around  14  K.  Solid  state  theories  usually  contained  in  textbooks  are  inadequate 
to  describe  the  physical  properties  (e.g.,  vibrational  excitations)  of  this  "quantum  solid."  About  20 
years  ago,  solid  hydrogen  was  widely  studied  precisely  because  the  large  zero-point  anharmonic 
vibrational  motion  of  the  lattice  molecules  imbued  this  system  with  interesting  properties.  These 
quantum  mechanical  effects  are  extremely  difficult  to  treat  from  the  point  of  view  of  dynamical 
computer  simulation.  Why  is  this?  Consider  the  case  of  the  impurity  self-diffusion.  The  diffusion 
constant  can  be  obtain  from  either  of  two  formulas: 

D  =  i  lim^(\qit)-q(0f}  (2) 

o  t-^oocLt 
or 

D  =  |j2r<q(r)-q(0)>  (3) 

where  the  notation  A{t)  s  denotes  a  quantum  Heisenberg  operator  for  any  operator 

"A".  Both  of  the  above  functions  are  many-body  quantum  dynamical  correlation  functions  [or,  in 


the  case  of  Equation  2,  can  simply  be  related  to  one].  Similarly,  the  intrinsic  dynamical 
recombination  step  of  two  impurities  can  be  described  by  "product  state"  population  correlation 
functions. 

In  a  quantum  solid  such  as  hydrogen,  the  above  quantum  time  correlation  functions  cannot  be 
approximated  using  the  classical  Newtonian  molecular  dynamics  (MD)  method.  The  centerpiece  of 
the  final  two  years  of  the  project,  therefore,  became  the  development  and  application  of  a  new 
method  developed  by  Cao  and  Voth^^-^o  for  computing  quantum  dynamical  time  correlation 
functions  such  as  those  above.  In  fact,  it  is  the  strength  of  this  method,  called  Centroid  Molecular 
Dynamics  (CMD),  which,  for  the  first  time,  made  it  possible  to  study  the  quantum  dynamics  of 
quantum  solids,  with  and  without  impurities.  This  method  is  now  making  it  possible  to  study  the 
self-diffusion  and  recombination  processes  of  impurities  in  cryogenic  HEDM  via  computer 
simulation  and,  thereby,  to  probe  the  microscopic  origins  of  HEDM  instabilities. 

This  research  project  has  been  both  critical  and  timely  for  the  Air  Force  HEDM  Program.  As 
the  HEDM  program  moves  into  the  future,  the  methodology  developed  during  the  Air  Force 
contract  can  be  expected  to  play  an  important  role  in  determining  the  feasibility  of  "scaling-up" 
certain  target  materials.  With  the  various  Department  of  Defense  (DoD)  High  Performance 
Computing  Centers  becoming  increasingly  functional,  the  methods  developed  with  the  support  of 
this  contract  are  being  implemented  on  extraordinarily  powerful  parallel  computing  platforms.  This 
should,  in  turn,  lead  to  condensed  phase  HEDM  simulations  of  an  unprecedented  scale  and  scope, 
while  addressing  some  very  challenging  problems  in  condensed  matter  theory. 

1.2  Brief  Summary  of  the  Research  Results 

•  The  most  important  result  obtained  during  the  contract  period  was  the  first  calculation  of  the 
quantum  mechanical  barrier  to  recombination  of  two  atomic  lithium  impurities  in  solid  hydrogen  at 
4  K.  The  barrier  was  found  to  be  on  the  order  of  130  K,  raising  the  theoretical  possibility  of 
trapping  lithium  impurities  in  solid  hydrogen  at  concentrations  of  over  1  mol  percent.  Though 
preliminary,  this  is  the  first  result  of  its  kind  and  is  encouraging  for  the  low  temperature  solid 
hydrogen  fuel  effort  of  the  Air  Force  HEDM  Program. 

•  A  computational  approach  to  study  low  temperature  HEDM  was  developed  called  CMD,  which 
incorporates  the  dominant  quantum  dynamical  effects  of  a  many-body  system  into  a  classical  MD 
framework.  A  powerful  CMD  algorithm  was  also  developed  called  Hyper-Parallel  CMD 
(HPCMD),  which  recently  achieved  3.1  billion  floating  point  operations  per  second  (GFLOPS) 
performance  with  80%  computational  efficiency  over  64  nodes  on  the  DoD  Maui  High 
Performance  Computing  Center  IBM  SP2.21 
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•  As  a  test  of  the  CMD  methodology,  the  quantum  self-diffusion  constant  for  liquid  para-hydrogen 
was  computed  for  several  temperatures  and  successfully  compared  by  experiment.  These 
simulations  demonstrated  both  the  accuracy  of  the  method  and  scalability  of  the  CMD  algoiithm.22 

•  The  quantum  phonon  spectmm  for  para-hydrogen  at  4  K  was  also  computed  and  successfully 
compared  by  experiment.  These  results  demonstrated  the  ability  to  reliably  study  low  temperature 
hydrogen  systems  of  interest  to  the  Air  Force  HEDM  Program  through  largescale  parallel  CMD 
simulations. 

•  Two  analytical  models  for  atomic  impurities  in  solid  hydrogen  were  developed,  which  have  been 
used  to  predict  the  equilibrium  structure  and  thermodynamic  properties  of  such  HEDMs.23,24 

•  The  nature  of  the  impurity  sites  for  lithium  atoms  in  solid  para-H2  and  ortho-D2  were  determined 
with  quantum  computer  simulation  methods.  The  electronic  excitation  spectrum  was  also  calculated 
and  found  to  agree  well  with  the  experimental  spectram  that  Fajardo  and  co-workers  obtained  at 
PL,  Edwards.25 

•  The  electron  spin  resonance  (ESR)  linewidth  for  a  hydrogen  atom  impurity  in  solid  p-H2  was 
calculated  from  first  principles  and  found  to  be  in  agreement  with  the  CMD  experiment.  The 
linewidth  was  predicted  to  be  the  same  for  substitutional  and  interstitial  impurity  sites  in  the  lattice. 
This  study  suggested  that  H-atom  impurities  in  solid  hydrogen  are  likely  to  be  quite  mobile  and, 
thus,  may  not  be  a  good  HEDM  candidate.26 

•  The  behavior  of  a  Li  atom  in  liquid  and  cluster  p-H2  was  studied  through  quantum  computer 
simulation  methods.^^ 

•  The  behavior  of  clusters  of  p-H2  were  studied  by  quantum  computer  simulation  methods.28-30 

2.0  RESEARCH  RESULTS,  DATA  OBTAINED,  AND  ACCOMPLISHMENTS 

2.1  Advances  in  Theoretical  Methodology 

Richard  Feynman  demonstrated  that  quantum  statistical  mechanics  could  be  reformulated  in 
terms  of  path  integrals. 3 1.32  For  example,  the  partition  function  for  a  system  is  given  by  the 
expression 
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Z  =  j  •  •  -  j  Dq(T)  exp{-S[q{x)y  fi)  (4) 

where  the  functional  integration  is  over  all  possible  paths  of  the  particles  q(T),  and  5[q(T)]  is  the 
imaginary  time  action  functional,  given  by 

•^[qW]  =  Jp^^/Tjiq^(T)-m-q(T)  +  V[q(T)]|  (5) 

While  the  numerical  evaluation  of  the  equilibrium  properties  of  a  condensed  phase  system 
using  Equation  4  might  seem  rather  daunting,  such  calculations  can  actually  be  readily 
performed.33-38  The  simplest  approach  involves  the  "primitive"  discretization  of  the  path  integral 
for  the  partition  function  into  P  "imaginary  time  slices,"  giving  a  formula  valid  for  sufficiently 
large  values  of  P,  i.e., 

N 

Z=  Urn  (6) 

j=i 

where  e  equals  P  /  P  and  the  isomorphic  effective  classical  potential  is  given  by 

p 

%(qb—qp)  =Z 

j=i 

and  the  subscript  "i"  denotes  the  coordinate  of  the  particle  at  imaginary  time  t  =  (t  -  ])he .  For  the 
calculation  of  properties  using  the  partition  function,  the  discretized  paths  are  cyclic  so  that  q, 
equals  qp+,.  The  primitive  discrete  representation  of  the  path  integral  is  isomorphic,  with  a  classical 
partition  function  for  a  polymer  composed  of  particles  located  at  the  positions  {q,}  (Fig.  1).  As 
can  be  deduced  from  Equations  6  and  7,  the  equilibrium  properties  of  condensed  matter  systems 
having  general  potentials  can  be  readily  evaluated  on  a  computer  by  Metropolis  Monte  Carlo  (MC) 
techiiiques.39  Many  studies  have  now  appeared  in  the  literature  in  which  path  integral  Monte  Carlo 
(PIMC)  or  path  integral  Molecular  Dynamics  (PIMD)  have  been  employed  to  simulate  a  quantum 
particle  (or  particles)  in  very  complex  environments  involving  hundreds  or  even  thousands  of 
atoms  and/or  molecules.33-38  Interestingly,  there  had  been  only  one  other  path  integral  simulation 
of  solid  hydrogen^o  before  those  carried  out  under  the  present  contract.  Moreover,  the  latter 
simulations  were  the  only  ones  to  study  impurities  in  solid  hydrogen.  Indeed,  it  is  largely  the 
strength  of  this  path  integral  simulation  approach  for  complicated  many-body  systems  that  allowed 
the  first  direct  calculations  of  a  Li  atom  trapped  in  a  solid  P-H2  matrix. 


mP 


2.2  (qf-qi+i)^+  ^^(qf) 


(7) 
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Figure  1 

A  Schematic  of  Two  Quantum  Particles  Interacting  Through  a 
Pairwise  Potential  in  the  Discretized  Path  Integral  Representation 


In  Figure  1,  each  quantum  particle  is  represented  by  an  isomorphic  classical-like  "polymer" 
of  P  quasiparticles  having  harmonic  intrapolymer  interactions  and  pairwise  interpolymer 
interactions.  The  value  of  the  discretization  parameter  P  in  this  case  is  6,  while  a  value  of  50-100  is 
typically  used  to  represent  each  quantum  p-PSQ,  molecule  in  the  spherical  (7=0)  approximation  at  a 
temperature  of  3-5  K. 

Whereas  the  research  supported  by  the  contract  over  the  first  three  years  was  focused  on  the 
structure  of  the  trapping  sites  for  impurities  in  solid  hydrogen,  research  in  the  last  two  years  began 
to  directly  address  the  dynamical  problem  of  impurity  diffusion  and  recombination  in  low 
temperature  HEDMs.  The  basic  components  of  this  problem  were  outlined  in  Section  1.1  (Aims  of 
Research).  If  such  dynamical  studies  were  to  be  carried  out,  however,  they  would  need  to  be 
quantum  dynamical  studies  [i.e.,  involving  the  effect  of  the  quantum  time  evolution  operator 
exp(— iTTr/S)].  For  low  temperature  HEDM  systems,  this  represented  an  insurmountable 
challenge,  since  to  solve  the  time-dependent  Schrddinger  equation  for  a  low-temperature,  nonlinear 
many-body  system  was  essentially  impossible  (at  least  exponentially  difficult!).  Fortunately,  a  key 
theoretical  breakthrough  allowed  study  through  computer  simulation  of  the  quantum  dynamics  of 
low  temperature  solids,  and  is  now  leading  to  direct  studies  of  impurity  diffusion  and 
recombination  in  cryogenic  HEDM.  This  breakthrough  is  called  Centroid  Molecular  Dynamics.^^- 
20 

In  "standard"  equilibrium  path  integral  simulations,^^'^*  the  isomorphic  quasiparticle 
polymer  (Fig.  1)  representing  a  quantum  particle  is  spread  out  in  space,  but  it  collapses  towards  a 
point  particle  and  becomes  more  classical-like  in  the  limit  of  high  temperature  or  large  mass.  It 
seems  reasonable,  therefore,  that  the  most  classical-like  variable  for  describing  the  equilibrium 
configurations  of  a  quantum  particle  must  be  the  path  centroid  variable,  '^1-47  defined  in  the 
discretized  path  integral  formalism  by 
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(8) 


1  4 

qo=i;2-q, 


PU 

where  are  the  coordinates  of  the  collection  of  P  quasiparticles  that  represent  the  quantized 
physical  particle  [Equations  6  and  7].  In  the  continuous  path  picture,  the  centroid  variable  is  given 
by 

1  1  ,*/? 

dr  q(T)  (9) 


]•  1  'V  ^ 


t=i 


From  the  latter  equation,  it  is  easy  to  see  that  as  >  0  (i.e.,  as  — >  0  or  T  — >  «>)  the  Feynman 
paths  that  represent  the  delocahzed  quantum  particle  will  shrink  down  to  a  classical  point  particle. 
However,  when  the  system  is  fully  quantized,  the  path  centroid  will  be  a  "ghost"  variable  within 
the  context  of  quantum  Boltzmann  statistical  mechanics. 

Before  describing  the  results  on  the  dynamical  implications  of  the  path  centroid  variable,  ^ '^■20 
it  is  important  to  note  that  the  interesting  equilibrium  properties  of  the  centroid  were  first 
discovered  and  exploited  by  Feynman  in  his  analysis  of  the  quantum  corrections  to  the  classical 
partition  function.4i-42  in  that  work,  Feynman  defined  a  path  integral  "centroid  density,"  given  by 

PMc)=  I  •••/  DqiT)5iqc  -  qo) expi-Slqir)]/ n}  (10) 


with  the  goal  of  deriving  a  simplified  quasiclassical  partition  function.  This  equilibrium  density  is 
calculated  by  fixing  the  centroid  variable  in  Equation  10  to  a  point  in  space  q^  and  then  integrating 

over-all  path  fluctuations  about  the  fixed  centroid.  In  the  discretized  path  integral  picture,  one 
instead  integrates  over-all  configurations  of  the  discretized  quasiparticles  having  a  fixed  centroid, 
subject  to  the  weighting  by  their  effective  Boltzmann-like  configurational  factor.  In  either  case 
(continuous  or  discretized),  the  centroids  of  the  Feynman  paths  are  located  at  the  position  in  space 
denoted  by  ^  "centroid  statistical  mechanics"  can  then  be  defined  based  on  the  centroid 

potential  ^^(qc),  defined  from  the  path  centroid  density  by 

exp[-/3v;.(qc)]  =  j  •••/  -  qQ)exp{-S[q{T)V%}  (11) 


where  a  normalization  factor  has  been  set  to  unity.  The  quantum  partition  function  Z  can  be 
determined  from  Equation  1 1  by  integrating  over-all  possible  centroid  positions  (i.e.,  the  centroid 
trace)  such  that 

Z  =  =  (12) 

Because  of  the  nature  of  the  path  centroid  variable,  the  centroid  density  in  Equations  10  and  1 1  is 
the  quantum  equivalent  to  the  classical  Boltzmann  density  for  the  particle's  equilibrium 
configurations. 

While  the  centroid  picture  in  equilibrium  statistical  mechanics  is  interesting  in  its  own  right,  it 
was  the  discovery  1 1 8  and  subsequent  derivation of  the  dynamical  properties  of  the  centroid 
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variable  that  have  made  research  on  impurity  diffusion  and  recombination  in  low  temperature  solids 
possible.  The  basic  result  of  CMD  is  that  time  correlation  functions  for  quantum  particles  in  many- 
body  systems  can  be  approximated  by  centroid  time  correlation  functions  obtained  from  running 
trajectories  on  the  effective  quantum  centroid  potential  (cf.  Fig.  2). 


In  Figure  2,  the  two  centroid  quasiparticles  interact  through  the  effective,  temperature- 
dependent  centroid  potential. 

The  centroid  trajectories  q^(0  in  CMD  are  generated  by  the  CMD  equation 

m-q,(t)  =  -  V,y,(q,)  (13) 

where  m  is  the  diagonal  particle  mass  matrix  and  the  effective  centroid  potential  is  defined 
according  to  Equation  1 1  as 

=  -  ^bTIApMc)] 

The  CMD  method  provides  a  way  to  directly  and  efficiently  include  quantum  zero-point  energy 
and  tunneling  in  molecular  dynamics  simulations.  It  has  also  been  showni^d9  that  the  self¬ 
diffusion  constant  can  be  calculated  from  CMD  by  studying  the  behavior  of  the  centroid  mean- 
squared  displacement  at  long  times  or,  equivalently,  by  a  Green-Kubo-like  relationship  for  the 
centroid  velocity  correlation  function  according  to  the  CMD  equivalents  of  Equations  2  and  3. 

One  remarkable  feature  of  CMD  is  that  it  allows  quantum  particle  motion  to  be  studied  with  a 
numerical  effort  that  scales  with  system  size  in  the  same  way  as  a  classical  MD  simulation.20 
Therefore,  the  dynamics  of  quantized  particles  in  general  many-body  systems  can  now  be 
simulated  without  having  to  solve  the  many-body  time-dependent  Schrodinger  equation.  A  number 
of  powerful  and  flexible  algorithms  for  solving  the  basic  CMD  equations  have  also  been 
developed.20  This  task  is  not  entirely  trivial,  since  the  centroid  potential  Equation  14  is  a 

kind  of  quantum  potential  of  mean  force  that  seemingly  requires  an  equilibrium  path  integral 
averaging  over  the  fixed-centroid  path  fluctuations  at  each  timestep.  A  "direct"  algorithm^o  for 
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integrating  the  CMD  equations,  which  performs  "on  the  fly"  the  statistical  averaging  implicit  in  the 
centroid  force  calculation,  combines  MD  with  MC  moves.  In  this  approach,  the  "natural"  CMD 
timestep,  being  on  par  with  the  timestep  in  the  classical  limit,  is  broken  into  Nf^Q  smaller 

timesteps.  At  each  of  these  centroid  configurations,  a  pass  of  pairwise  MC  moves  of  the 
quasiparticles  is  carried  out  to  sample  the  discretized  path  configurations  while  enforcing  the 
centroid  constraint.  The  centroid  forces  are  then  computed  and  the  centroids  moved  according  to 
those  forces  within  the  velocity  Verlet  algorithm  for  the  small  timestep.  This  algorithm  has  proven 
to  be  efficient  for  numerically  solving  the  CMD  equations. 

2.2  Path  Integral  Simulations  of  Impurity  Trapping  and  Spectroscopy  in  Quantum 
Solids,  Liquids,  and  Clusters 

Experiments  at  Phillips  Laboratory,  Edwards  AFB,  have  demonstrated  that  a  low 
concentration  of  lithium  atoms  can  be  metastably  trapped  in  both  solid  H2  and  D2  in  the 
temperature  range  r<  5  K.b2  This  is  of  interest  because  the  ability  of  hydrogen  to  act  as  a  rocket 
fuel  might  be  significantly  enhanced  if  small  quantities  of  such  light  impurities  are  introduced  into 
the  system.  To  better  understand  the  properties  of  such  metastable  quantum  "alloys,"  it  was  a 
priority  to  characterize  the  trapping  sites  of  individual  Hthium  atoms  in  the  solid  hydrogen  matrix. 

As  a  first  step,  PIMC  studies  of  solid,23.24  liquid,^^  and  cluster^^.so  para-hydrogen  and 
ortho-deuterium  were  performed  using  a  constant  pressure  version  of  the  path  integral  formulation 
described  in  the  previous  section.  A  semi-empirical  intermolecular  pair  potential,  which  had  been 
obtained  by  others  a  number  of  years  ago  through  a  fitting  to  solid  state  data,^^  was  used  in  our 
studies  for  the  hydrogen.  This  potential  gives  a  rather  good  account  of  pV  compression  data.  [It 
should  be  noted  that  a  pressure  of  only  200  MPa  (2  kbar)  is  sufficient  to  change  the  molar  volume 
of  solid  hydrogen  from  23.2  cc/mol  to  ca.  17  cc/moL]  The  calculated  energy  and  molar  volume  of 
the  hydrogen  liquid  from  the  PIMC  simulations  were  found  to  be  in  excellent  agreement  with 
experiments  for  temperatures  ranging  from  the  triple  point  up  to  the  boiling  point.  These  results 
justified  the  choice  of  the  intermolecular  potential  used  in  the  simulations,  as  well  as  the 
computational  methodology. 

Once  the  reliability  of  the  PIMC  code  had  been  established  for  the  case  of  solid  and  liquid 
hydrogen,  the  simulation  effort  focused  on  a  lithium  atom  impurity  immersed  in  bulk  liquid^^  and 
solidus  para-hydrogen,  as  well  as  ortho-deuterium.  A  question  first  arose  as  to  the  potential  for  the 
lithium  (guest)  -  hydrogen  (host)  interaction.  A  lithium-hydrogen  pair  potential  was  constmcted  by 
fitting  the  results  of  ab  initio  quantum  chemical  c2ilculations  and  is  shown  in  Figure  3.  This 
procedure  has  since  been  further  refined. 
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Figure  3 

A  Comparison  of  the  Hj— H2  and  Li-H2  Empirical  Pair  Potentials  Employed 
in  the  PIMC  Simulations  of  a  Li  Impurity  in  a  Solid  P-H2  Host 


Since  an  isolated  lithium  atom  is  much  larger  than  the  host  hydrogen  molecules,  trapping 
sites  consisting  of  one  to  six  vacancies  in  the  hydrogen  lattice  were  investigated  (see  Fig.  4). 
Interestingly,  all  of  the  sites  were  found  to  be  comparable  in  energy.  This  is  due  to  the  large 
compressibility  of  para-hydrogen  and  ortho-deuterium  solids,  which  permits  the  lattice  to  relax  to 
comfortably  accommodate  the  impurity. 


Figure  4 

A  Schematic  Showing  Various  H2  Molecules  Which  can  be  Removed 
from  hep  p-Hi  Lattice  to  Accommodate  an  Atomic  Li  Impurity 


The  numbering  in  Figure  4  refers  to  the  number  of  hydrogen  vacancies  and  their  position. 
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The  treatment  of  any  system  consisting  of  an  impurity  plus  a  solid  with  vacancies  requires 
some  care.  The  presence  of  the  vacancies  can  give  rise  to  a  number  of  different  local  minima  in  the 
solid  that  can  trap  the  impurity.  These  minima  can  be  separated  by  large  barriers  that  are  difficult  to 
overcome  on  the  "time  scale"  of  a  Monte  Carlo  simulation.  Therefore,  a  relatively  simple  and 
unbiased  PIMC  method  was  developed  for  the  determination  of  the  likely  configurations  of 
vacancies  bound  to  the  lithium  impurity  in  the  para-hydrogen  or  ortho-deuterium  systems.^s  In  one 
of  the  earlier  studies, a  lithium  impurity  in  liquid  para-hydrogen  at  T  =  14-20  K  was  investigated 
(Fig.  5).  For  such  a  liquid,  the  relevant  configurations  could  be  sampled  efficiently,  and  the 
solvation  shell  of  the  lithium  impurity  was  found  to  consist  of  approximately  20  molecules  (Fig. 
6).  Several  configurations  of  the  "solvated"  lithium  atom  and  its  first  20  neighbors  in  the  liquid 
were  then  saved  for  subsequent  use.  We  refer  to  such  configurations  as  "liquid  cores."  A  trapping 
site  in  the  hydrogen  solid  was  then  prepared  using  the  liquid  cores  in  the  following  manner.  An 
arbitrary  point  in  an  equilibrated  neat  solid  configuration  was  chosen  and  20+nv  molecules 
removed.  The  liquid  core  was  then  inserted  into  the  solid  and  a  constant  pressure  PIMC  simulation 
commenced.  Equilibration  to  a  relaxed  impurity  "implantation"  configuration  was  found  to  be  rapid 
(~  5000  steps).  This  procedure  was  repeated  for  a  number  of  different  liquid  cores  and  neat  solid 
configurations.  The  assumption  in  this  procedure  is  the  physically  reasonable  one;  i.e.,  that  the 
vacancies  occur  only  in  the  first  neighbor  shell  around  the  lithium  atom.  Interestingly,  for  a  given 
riv,  relatively  few  relaxed  configurations  were  observed.  Shown  in  Figures  7  and  8  are  the 
resulting  sites  with  =  3  and  4  vacancies  (cf.  Fig.  4)  that  were  energetically  favored.  These 
trapping  sites  always  contain  Li  in  asymmetric  locations,  i.e.  shifted  off  the  ideal  interstitial  or 
lattice  sites. 

The  inhomogeneously  broadened  dipole  spectra  of  the  lithium  impurity  in  the  various  sites 
were  next  calculated^^  and  compared  to  the  experimental  results  of  Fajardo. ^2  This  was 
accomplished  by  taking  equilibrium  configurations  from  the  constant  pressure  PIMC  runs, 
replacing  the  empirical  lithium  atom  with  an  electron  and  a  lithium  ion  core,  both  of  which  interact 
with  the  hydrogen  host  through  pseudopotentials,  and  then  determining  the  electronic  states  using  a 
radial  fast  Fourier  Transform  Lanczos  method.27  The  spectra  for  Li  in  various  sites  of  the 
hydrogen  solid  are  shown  in  Figure  9.  This  approach  accounts  for  the  zero-point  motion  and 
thermal  effects  present  in  the  system.  Based  on  these  calculations,  lithium  atoms  were  suggested  to 
preferentially  occupy  a  three-vacancy  trapping  site  in  para-hydrogen,  while  in  ortho-deuterium  a 
four-vacancy  trapping  site  seems  to  be  favored.  Interestingly,  the  high  energy  tail  of  the 
experimental  spectrum  is  underestimated  in  all  of  the  calculations.  The  "average"  energy  of  this  tail 
region  suggests  that  it  may  in  fact  be  caused  by  a  rotational  (/  =  0  — >  2)  Franck-Condon  transition 
in  the  neighboring  hydrogen  molecules;  such  effects  have  been  neglected  in  the  simulation  but 
could  be  included  in  the  future  through  an  explicit  treatment  of  the  hydrogen  rotations. 
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Figure  5 

A  Snapshot  Taken  from  a  Constant  Pressure  PIMC  Calculation 
of  a  Lithium  Atom  in  Liquid  Hydrogen 


Each  blue  sphere  in  Figure  5  represents  a  para-hydrogen  molecule  and  the  larger  yellow  ball 
the  lithium  impurity.  The  closest  20  neighbors  to  the  impurity  are  colored  light  blue.  The  latter 
constitute  the  “liquid  core”  that  is  “implanted”  into  hep  solid  hydrogen  as  a  means  of  searching  for 
the  possible  trapping  sites  of  ablated  lithium  atoms. 


Figure  6 

The  Solvation  Shell  of  the  Lithium  Impurity  in  Liquid  Para- 
Hydrogen  Consisting  of  Approximately  20  Hydrogen  Molecules 
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Figure  7 

Top  View  of  a  Configuration  Taken  from  a  Constant  Pressure  PIMC 
Calculation  of  a  Lithium  Atom  Trapping  Site  in  Solid  Hydrogen 


The  impurity  site  in  Figure  7  was  created  by  first  removing  23  host-lattice  hydrogen 
molecules  and  then  “implanting”  a  “liquid  core“  consisting  of  a  single  hthium  atom  and  20 
hydrogen  molecules.  The  asymmetric  location  of  the  lithium  atom  in  an  otherwise  perfect 
hexagonal  lattice  with  three  vacant  sites  was  the  result  of  the  constant  pressure  PIMC  calculation. 
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Figure  8 

Top  View  of  a  Configuration  Taken  from  a  Constant  Pressure  PIMC 
Calculation  of  a  Lithium  Atom  Trapping  Site  in  Solid  Hydrogen 

The  impurity  site  in  Figure  8  was  created  by  first  removing  24  host-lattice  hydrogen 
molecules  and  then  "implanting"  a  "liquid  core"  consisting  of  a  single  lithium  atom  and  20 
hydrogen  molecules.  The  asymmetric  location  of  the  lithium  atom  in  an  otherwise  perfect 
hexagonal  lattice  with  four  vacant  sites  was  the  result  of  the  constant  pressure  PIMC  calculation. 
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Figure  9 

The  Calculated  Dipole  Spectrum  of  a  Lithium  Atom  in  Solid  Para- 
Hydrogen  at  T  =  4  K  for  Various  Trapping  Sites 

The  trapping  sites  (cf.  Figs.  3,  6  and  7)  are  compared  above  to  the  experimental  data  of  Fajardo.l’^ 
Complementary  calculations  on  the  Li/H2  system  were  performed  using  the  variational 
quantum  Einstein  modeP^  for  selected  trapping  sites  (see  Section  2.3,  Theoretical  Description  of 
Impurities  in  Solid  Hydrogen).  Based  on  free  energy  considerations,  it  was  found  that  the  Einstein 
model  favors  the  four-vacancy  site  in  both  para-hydrogen  and  in  ortho-deuterium.  Since  the 
Einstein  model  is  a  local  harmonic  one,  its  difference  from  the  simulation  results  illustrates  the 
importance  of  anharmonic,  collective  quantum  motion  in  determining  the  stracture  of  the  impurity 
trapping  site. 

The  interaction  of  Li  with  clusters  containing  either  12,  13,  32,  33,  or  34  para-hydrogen 
molecules  were  also  investigated.^^  The  lithium  atom  was  found  to  reside  on  or  near  the  cluster 
surface  even  though  the  clusters  were  studied  at  much  lower  temperatures  than  the  liquid;  i.e.,  2.5 
to  6  K.  Figure  10  shows  representative  configurations  taken  from  the  path  integral  calculations  for 
the  cluster  with  33  para-hydrogen  molecules.  Although  in  ail  cases  the  Li  atom  resides  on  the 
outside  of  the  cluster,  perturbations  of  the  stracture  are  observed  in  comparison  to  neat  para- 
hydrogen  clusters,  which  were  also  studied.^^  For  both  the  clusters  and  the  neat  liquid,  the 
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inhomogeneously  broadened  dipole  spectrum  of  the  lithium  atom  was  again  calculated  using  the 
radial  fast  Fourier  Transform  Lanczos  method.  In  the  clusters,  the  spectra  exhibit  a  main 
absorption  band  near  the  unperturbed  atomic  Li  value  and  a  second  asymmetric  band  shifted  to  the 
blue.  The  latter  is  identified  as  arising  from  the  p-orbital  oriented  radially  towards  the  cluster, 
whereas  the  main  band  is  composed  of  the  two  p-orbitals  oriented  parallel  to  the  cluster  surface. 
The  spectrum  of  Li  in  the  liquid  is  significantly  broader  than  the  cluster  spectra.  The  ionization 
spectrum  of  Li  attached  to  the  para-hydrogen  clusters  was  also  smdied.  The  spectra  progressively 
red-shift  and  broaden  from  the  atomic  value  with  increasing  cluster  size. 


Figure  10 

A  Snapshot  from  a  PIMC  Simulation  of  the  Li(p-H2)33  Cluster 


In  Figure  10,  note  the  ordering  of  the  hydrogen  molecules  and  the  fact  that  the  lithium  atom  is 
on  the  surface  of  the  cluster. 

2.3  Theoretical  Description  of  Impurities  in  Solid  Hydrogen 

Another  research  effort  involved  the  development  of  a  theoretical  framework  by  which  to 
characterize  and  understand  the  structure  and  thermodynamics  of  impurities  in  the  P-H2  solid.23.24 
These  analytical  developments  were  variational  in  nature  and  based  on  the  well-known  Gibbs- 
Bogoliubov  (GB)  inequality.  This  inequality  states  that  the  true  Helmholtz  free  energy  of  a  system 
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is  the  lower  bound  of  the  sum  of  the  free  energy  for  a  variational  "reference  system"  and  the 
difference  between  the  total  energies  of  the  exact  and  reference  systems,  averaged  over  the 
canonical  ensemble  generated  from  the  reference  system.  Explicitly,  the  GB  inequality  is  given  by 

F  <  FiO  =  +  {H-  H,,f{0)ref  (15) 

where  C  denotes  the  set  of  variational  parameters.  Variational  theories  based  on  the  above  equation 
involve  the  specification  of  a  suitable  (and  tractable!)  reference  Hamiltonian,  which  depends  on  one 
or  more  variational  parameters,  and  the  subsequent  minimization  of  the  right-hand  side  of  Equation 
15  with  respect  to  those  parameters. 

As  a  first  step  in  the  theoretical  analysis,  a  variational  Einstein  model  for  describing  low 
temperature  solids,  with  and  without  impurities,  was  developed  from  the  Feynman  path  integral 
perspective. 23  The  theory  can  be  used  to  predict  fully  quantum  mechanical  values  for  the 
thermodynamics  (e.g.,  free  energy,  entropy,  internal  energy,  etc.)  and  the  equilibrium  structure 
(e.g.,  pair  and  angular  correlation  functions)  of  a  solid.  The  independent  harmonic  oscillator 
assumption  implicit  in  the  Einstein  model  allows  the  results  to  be  cast  in  a  straightforward  analytic 
form.  Additionally,  the  path  integral  formulation  of  the  model  yields  solutions  which  explicitly 
depend  on  the  path  integral  discretization  parameter  P.  One  can  thus  systematically  examine  the 
equilibrium  behavior  of  a  solid,  ranging  from  the  classical  to  the  quantum  limits.  The  Einstein 
model  was  applied  to  examine  the  behavior  of  solid  hydrogen  and  solid  hydrogen  containing 
lithium  impurities. 

Although  the  description  of  solid  hydrogen  with  the  path  integral  Einstein  model  agreed  fairly 
well  with  simulations  in  the  nearly  classical  limit,  somewhat  larger  discrepancies  were  observed  in 
the  fully  quantum  limit.23  This  difficulty  is  at  least  partially  a  result  of  the  independent  harmonic 
oscillator  assumption  for  the  solid.  Since  the  excitations  in  a  quantum  solid  are  mainly  acoustic 
phonons,  the  particles'  motions  are  expected  to  be  strongly  correlated  with  one  another,  and  the 
independent  harmonic  oscillator  assumption  may  become  inappropriate  to  characterize  solids  at  low 
temperatures.  The  acoustic  phonons  can,  however,  be  correctly  characterized  by  a  Debye  model, 
the  parameters  of  which  could  also  be  determined  variationally  from  Equation  15. 

Impurities,  on  the  other  hand,  destroy  the  translational  symmetry  of  a  solid,  so  the  dynamical 
picture  for  a  solid  with  impurities  becomes  greatly  complicated.  The  motion  of  an  impurity  can 
either  be  localized  when  its  vibrational  frequencies  are  much  higher  than  the  Debye  frequency,  or 
the  impurities'  motion  may  be  strongly  mixed  with  the  delocalized  vibrational  states  if  its 
vibrational  frequencies  are  smaller  than  the  Debye  frequency.  In  general,  the  motions  of  impurities 
are  represented  by  neither  of  those  two  extreme  cases.  The  theoretical  effort,  therefore,  turned  to 
the  development  of  an  analytical  model  that  bridges  the  two  extremes  and  accurately  accounts  for 
the  effects  of  impurities  in  a  low  temperature  solid.24  This  new  variational  "mixture"  model  was 
formulated  to  study  the  effect  of  substitutional  impurities  in  solids  at  very  low  temperature.  Similar 
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to  the  Einstein  model,  the  mixture  model  is  a  self-consistent  harmonic  one.  In  the  latter  model, 
however,  except  for  the  impurities,  the  deviations  of  all  the  particles  from  their  thermodynamical 
equilibrium  positions  are  assumed  to  arise  from  the  vibrations  of  an  effective  harmonic  system  in 
which  each  particle  is  coupled  to  all  others.  To  account  for  the  nature  of  the  motion  of  impurities, 
additional  degrees  of  freedom,  which  are  used  for  the  localized  modes,  were  introduced  into  the 
model.  The  application  of  the  GB  variational  principle  was  carried  out  in  the  extended  space  in 
order  to  determine  the  variational  parameters.  The  actual  motions  of  impurities  were  then 
partitioned  between  the  correlated  motion  with  the  lattice  and  the  independent  harmonic  oscillations 
via  a  projection  scheme.  This  theoretical  analysis,  which  was  not  straightforward,  effectively 
"mixes"  the  Einstein  and  Debye  models  variationally,  allowing  for  extended  lattice  vibrations  and 
arbitrarily  localized  impurity  modes.  The  averaged  effect  of  the  anharmonicity  in  the  pair  potentials 
is  included  self-consistently  in  the  variational  model,  and  the  impurities  are  allowed  to  induce 
structural  distortions  at  constant  pressure.  The  (quite  accurate)  predictions  of  the  mixture  model  for 
the  pair  distribution  function  of  pure  P-H2,  and  for  H-H2  the  distribution  for  solid  hydrogen  with  a 
substitutional  atomic  hydrogen  impurity,  are  shown  in  Figure  1 1  and  compared  there  to  exact 
PEMC  results. 


r  (a.u.) 


Figure  11 

Pair  Distribution  Functions  in  Solid  P-H2 


The  top  panel  of  Figure  1 1  shows  the  theoretical  "mixture  model"  prediction  of  the  H2-H2 
pair  distribution  function  (solid  line)  and  the  H-H2  pair  distribution  function  (dashed  line)  in  a  p- 
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H2  solid  with  a  substitutional  atomic  hydrogen  impurity.  The  lower  panel  shows  the  same 
quantities  as  calculated  from  a  PIMC  simulation. 

2.4  Structure  and  Spectroscopy  of  Hydrogen  Atoms  in  Solid  Molecular  Hydrogen 

For  a  number  of  years,  there  has  been  considerable  interest  in  the  structure  and  dynamics  of 
atomic  hydrogen  impurities  in  solid  hydrogen^^'^^  because  of  their  theoretical  potential  to 
significantly  enhance  the  characteristics  of  rocket  fuels.  For  example,  the  ESR  absorption  spectra 
of  the  trapped  H  atoms  in  solid  hydrogen  has  been  reported.^^  These  experimental  results 
presented  a  good  challenge  for  a  theoretical/computational  study  The  main  purpose  of  this  effort 
was  to  derive  a  formula  for  the  ESR  linewidth  that  takes  into  account  the  localized  nature  of  H  atom 
impurity  electron  (Fig.  12)  and  the  large  fluctuations  of  the  quantized  nuclear  configurations.  The 
modified  formula  would  then  be  employed  in  a  first-principles  calculation  to  study  the  effects  of  the 
local  distortions,  as  well  as  the  molecular  hydrogen  vibrations  and  rotations,  on  the  linewidth  of 
the  ESR  spectra. 


Figure  12 

A  Schematic  Sketch  of  the  Square  of  the  Localized  Atomic  Wave  Function, 

4>(r),  of  an  H  Atom  and  a  One-Electron  Wave  Function  of  an  H2  Molecule 

In  Figure  12,  the  solid  line  represents  the  atomic  wave  function  of  the  H  atom,  the  short 
dashed  line  corresponds  to  an  H2  orientation  parallel  to  the  line  which  connects  the  H  atom  and  the 
center  of  mass  of  the  H2  molecule,  and  the  long  dashed  lines  is  for  an  H2  orientation  perpendicular 
to  the  line. 
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The  first  issue  considered  was  the  nature  of  the  electronic  wave  function  of  the  trapped 
impurity.  When  an  H  atom  is  trapped  in  solid  hydrogen,  according  to  the  Pauli  principle,  its  atomic 
wave  function  must  adjust  itself  to  avoid  overlap  with  the  wave  functions  of  the  molecular  H2 
electrons  that  have  the  same  spin  (the  molecular  orbitals  of  the  solid  molecules  must  change  as 
well).  Based  on  the  requirement  that  the  total  electronic  energy  must  be  a  minimum,  an  equation 
was  derived  from  Hartree's  self-consistent  field  theory  to  uniquely  determine  the  orthogonalization 
of  atomic  orbitals  (Fig.  13).  Then,  based  on  the  perspective  that  hyperfine  interactions  between  the 
electron  of  the  impurity  and  the  nuclear  moments  of  adjacent  particles  lead  to  the  ESR  hnewidth,  a 
formula  was  derived  that  includes  the  ensemble  average  over  the  quantized  nuclear 

configurations.^^ 


Figure  13 

The  Correction  of  the  Atomic  Wave  Function,  4>(r),  for  an  H 
Atom  Caused  by  the  Presence  of  an  H2  Molecule 

In  Figure  13,  the  H2  molecule  is  R  =  6.5  atomic  units  (a.u.)  away  with  its  orientation 
pointing  toward  the  H  atom.  The  solid  line  corresponds  to  our  rigorous  orthogonalization 
procedure,  while  the  long  dashed  line  is  the  atomic  wave  function  obtained  when  an  approximate 
Gram-Schmidt  orthogonalization  is  used.  The  short  dashed  line  is  the  original  zeroth-order  atomic 
wave  function.  Note  the  anti-bonding-like  behavior  of  the  orthogonalized  hydrogen  wavefunction. 

The  newly  derived  formula  was  then  employed  in  the  computation  of  the  linewidths  of  the 
ESR  spectra  for  impurities  trapped  in  substitutional  and  interstitial  sites  in  solid  hydrogen.  The 
nuclear  configurations  were  generated  from  constant  pressure  PIMC  simulations.  The  numerical 
calculations  showed  that  the  inclusion  of  zero-point  vibrations  drastically  increased  the  linewidth  of 
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the  ESR  spectra  over  the  linewidth  calculated  simply  from  the  average  equilibrium  configuration. 
From  the  simulations,  it  was  also  found  that  an  H  atom  in  an  interstitial  site  will  effectively  push  its 
nearest  neighbors  outward.  More  importantly,  when  the  quantum  effects  of  the  solid  are  fully  taken 
into  account,  the  lattice  is  not  able  to  maintain  the  cavity  of  the  interstitial  impurity,  and  the  local 
structure  around  the  impurity  becomes  the  same  as  that  for  a  substitutional  impurity  (Fig.  14). 
Because  of  this  large  structural  relaxation,  the  linewidth  for  the  H  atom  trapped  at  interstitial  sites 
decreases  to  the  same  value  as  that  for  an  H  atom  trapped  in  substitutional  sites.  The  calculated 
ESR  linewidth  in  either  site  was  found  to  be  in  virtually  exact  agreement  with  the  experimental 
result  of  0.34  G.  This  conclusion  remains  unchanged  when  a  pressure  as  high  as  10  kbar  is 
applied  to  the  solid.  Since  the  observed  lattice  relaxation  and  calculated  hnewidths  suggest  that  all 
trapping  sites  become  essentially  equivalent  (i.e.,  somewhat  like  a  fluid),  this  theoretical  study  may 
cast  some  doubt  on  the  ability  to  trap  and  stabilize  a  significant  mol  fraction  of  atomic  hydrogen 
impurity  in  solid  p-H2. 


Figure  14 

A  Comparison  of  the  Pair  Correlation  Function,  g(r).  Between 
an  H  Atom  and  P-JI2  Molecules 

The  curve  in  Figure  14  was  generated  from  a  PIMC  calculation  at  5  K  with  P  =  50  and  for 
different  initial  conditions.  The  solid  line  corresponds  to  the  H  atom  initially  in  an  interstitial  site 
while  the  dashed  line  is  for  it  in  a  substitutional  site. 
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2.5  Centroid  Molecular  Dynamics  Algorithms 

To  carry  out  the  quantum  dynamical  CMD  studies  of  low  temperature  HEDM  systems, 
efficient  and  effective  algorithms  were  needed These  are  described  in  the  next  paragraph. 

The  discretized  path  integral  in  Equation  6  converges  to  the  true  result  as  P  — >  oo.  In  practice, 
a  finite  value  of  P  is  chosen  to  ensure  that  the  quantities  of  interest  are  adequately  converged.  For 
example,  solid  para-hydrogen  has  pronounced  quantum  effects  at  the  temperatures  studied,  and  a 
large  value  of  P  is  necessary  at  the  lowest  temperatures  (i.e.,  P  =  50-100  at  5  K)  to  obtain 
convergence  of  the  dynamical  and  structural  properties.  Thus,  a  reasonable  amount  of  centroid- 
constrained  averaging  is  necessary  at  each  CMD  timestep  to  calculate  the  centroid  force,  making  its 
calculation  computationally  challenging.  However,  this  averaging  lends  itself  to  quick  and  efficient 
evaluation  on  a  parallel  computer. 

Normally,  for  a  direct  CMD  calculation^®  carried  out  in  serial  fashion,  one  PIMD  (or, 
alternatively,  PIMC)  trajectory  is  run  for  the  number  of  time  steps,  M,  to  evaluate  the  centroid 
force  at  each  CMD  timestep  (M  is  chosen  based  on  convergence  issues).^®  In  the  HPCMD 
implementation,^!  centroid  force  is  calculated  in  parallel  by  averaging  m  independent  PIMD 
trajectories,  and  each  trajectory  is  run  for  M/m  timesteps  (Fig.  15).  In  addition  to  being  faster,  this 
procedure  in  principle  samples  phase  space  more  efficiently  than  the  traditional  method  of 
evaluating  the  centroid  force  with  a  single  PIMD  simulation  because  the  configurations  are  not  as 
highly  correlated.  Each  PIMD  centroid  averaging  trajectory  is  started  from  the  same  initial 
quasiparticle  configuration,  but  the  velocities  are  randomized  at  the  beginning  of  the  simulation.  A 
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Figure  15 

A  Schematic  Flow  Diagram  of  the  Parallelism  Inherent  in  the 
Hyper-Parallel  CMD  (HPCMD)  Algorithm 
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Nose  chain  can  be  linked  to  the  centroid-constrained  quasiparticles  in  PIMD  to  ensure  ergodic 
sanipling58>59  (chains  of  length  two  have  been  used  for  each  hydrogen  molecule).  Alternatively,  a 
staging  PMC  algorithm  can  be  used  for  the  centroid  force  averaging,  as  originally  described  in 
Reference  20,  which  can  lead  to  more  rapid  convergence  of  the  force.58,59 

The  method  of  running  independent  trajectories  in  parallel  to  obtain  the  centroid  force  is 
termed  the  "first  tier"  of  parallelism  in  HPCMD.  This  part  of  the  algorithm  can  be  run  over  any 
number  of  processors  on  a  parallel  computer,  up  to  the  point  where  one  PMD  timestep  is  run  per 
processor  node,  which  is  the  maximum  division  of  labor  over  M  processors,  at  this  level.  As 
outlined  in  Reference  21,  the  number  of  timesteps  M  can  be  adjusted  to  test  convergence  of  the 
centroid  force  and  the  uniqueness  of  the  centroid  trajectories. 

As  a  result  of  this  parallelism,  each  PMD  simulation  does  M/m  PMD  timesteps  to  evaluate 
the  centroid  force,  but  each  one  of  these  timesteps  can  also  be  done  in  a  parallel  manner  over 
imaginary  time  slices  with  minimal  inter-node  communication  and  perfect  load  balancing  (cf.  Fig. 
15).  Since  each  trajectory  has  N  physical  particles,  and  each  of  these  is  discretized  into  P 
equivalent  quasiparticles,  the  overall  force  calculation  scales  approximately  asN^xP  (or  less, 
depending  upon  the  range  of  the  interactions).  Thus,  the  obvious  parallel  solution  is  to  evaluate  the 
intermolecular  forces  (often  99%  of  the  computational  effort)  for  each  imaginary  time  slice  in 
parallel  at  each  PMD  timestep.  This  parallel  evaluation  of  the  intermolecular  force  at  each  PMD 
step  is  termed  the  "second  tier"  of  hyper-parallelism.  This  tier  of  the  algorithm  can  be  run  on  any 
number  of  processors  up  to  a  total  of  P,  the  number  of  imaginary  time  slices.  It  is  most  efficient 
when  P  is  an  integer  multiple  of  the  number  of  processors. 

With  the  two  levels  of  parallelism  just  outlined,  the  CMD  calculation  can  be  run  over  Mx  P 
processors.  Therefore,  the  maximum  division  of  labor  is  when  all  M  PMD  calculations  of  the 
centroid  force  at  each  CMD  timestep  are  done  on  separate  processors,  and  in  turn  all  P  imaginary 
time  slices  for  each  PMD  calculation  are  done  on  P  separate  processors.  In  practice,  however, 
this  degree  of  parallelism  may  not  always  be  utilized  because  of  the  limited  number  of  processors 
available. 

A  "third  tier"  of  hyper-parallelism  is  even  possible  for  very  large  systems  (cf.  Fig.  15).  The 
parallel  evaluation  of  the  intermolecular  force  for  each  PMD  step  can  still  be  an  problem  at  each 
imaginary  time,  slice  for  longer-ranged  forces.  This  force  evaluation  can  be  done  in  parallel,  but 
because  of  inter-node  communication  bottlenecks,  this  is  only  efficient  for  large  systems  (i.e., 
usually  1000  particles  or  more).  This  force  loop  parallelism  has  been  implemented  in  other 
studies^®’^!  but  not  yet  in  the  para-hydrogen  system,  as  the  intermolecular  force  problem  is  not 
exceptionally  expensive  until  the  number  of  particles  reaches  many  thousand. 

Benchmarks  of  the  HPCMD  algorithm  were  carried  out  for  liquid  para-hydrogen  on  an  IBM 
SP2.  The  performance  and  scalability  of  the  algorithm  is  shown  in  Figure  16.  The  plot  of  speedup 


versus  number  of  processors  on  the  SP2  was  almost  linear  over  a  range  from  1  to  64  nodes.  On  an 
IBM  590  (one  wide  node  of  the  SP2),  the  serial  version  of  the  code  performed  at  70  million 
floating  point  operations  per  second  (MFLOPS),  while  on  the  Cray  C90  the  serial  version  of  the 
code  ran  at  280  MFLOPS.  In  contrast,  the  parallel  code  ran  at  approximately  3.1  GFLOPS  over  64 
nodes.  At  this  speed,  the  code  utilized  the  parallel  CPU  array  with  about  an  80%  efficiency  (real 
speedup/maximum  theoretical  speedup),  and  it  ran  more  than  10  times  as  fast  as  the  single 
processor  on  a  Cray  C90.  The  simulations  given  in  Figure  16  were  carried  out  for  1440  particles, 
and  the  efficiency  would  clearly  increase  with  larger  system  sizes.  It  is  anticipated  that  for  very 
large  hydrogen  systems,  a  HPCMD  code  could  be  developed  that  achieves  50-100  GFLOPS 
performance  on  a  current-generation  IBM  SP2. 


Figure  16 

The  Scaling  of  the  HPCMD  Algorithm  with  Number  of  Parallel 
Nodes  on  an  IBM  SP2  Computer 

The  test  simulation  in  Figure  16  was  for  liquid  p-H2  at  14  K  with  1440  molecules  in  the  simulation 
cell. 

2.6  Quantum  Dynamical  Simulations  of  Liquid  Hydrogen 

As  an  early  test  of  the  HPCMD  algorithm,  velocity  auto-correlation  functions,  for 

liquid  para-hydrogen  were  computed.2f22  j^e  pair  potential  for  the  interparticle  interactions  was 
taken  to  be  that  of  Silvera  and  Goldman.'^^  The  path  integral  was  discretized  into  P  =  16 
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quasiparticles  at  25  K,  while  at  14  K  it  was  found  that  P=  50  was  necessary  for  the  structural 
properties  to  be  well  converged  24.25  jn  the  correlation  function  calculations,  it  was  found  that  180 
periodically  replicated  hydrogen  molecules  were  adequate  for  convergence  with  system  size. 

At  25  K,  independent  trajectories  with  lengths  of  4,  6,  8,  and  20  ps  were  used  to  test  the 
convergence  of  Cy^it).  There  were  only  minor  differences  between  the  correlation  function 

computed  from  an  average  of  two  8  ps  trajectories  and  that  computed  from  a  20  ps  trajectory.  The 
velocity  correlation  function  is  shown  in  Figure  17.  The  converged  CMD  self-diffusion  constant 
was  calculated  from  the  average  of  two  8  ps  trajectories.  The  CMD  method  reproduced  the 
experimental  self-diffusion  constant  quite  well  (1.54  ±  0.1  A^/ps  versus  1.6  A^/ps 
experimentally).  The  classical  MD  result  in  which  a  constant  pressure  classical  MD  simulation  was 
first  used  to  equilibrate  the  system  volume24.25  given  by  0.5  ±  0.05  A^/ps,  so  the  quantum 
effects  are  significant  even  at  this  temperature. 


Figure  17 

Plot  of  the  CMD  Velocity  Auto-Correlation  Functions  for  Liquid 
P-H2  at  25  K  and  V  =  31.7  cm®  mol  * 

The  quantum  CMD  value  of  the  self-diffusion  constant  is  found  to  be  1.54  A^  ps-^,  while  the 
experimental  value  is  1.6  ps'^  Calculations  at  other  densities  give  similar  agreement  between 
CMD  and  experiment. 

At  14  K,  two  9.8  ps  HPCMD  trajectories  were  run.  A  serial  calculation  of  these  trajectories 
by  direct  CMD  was  problematical  for  a  reasonable  turnaround  time  because  of  the  number  of 
quasiparticles  and  the  amount  of  averaging  required  to  calculate  the  centroid  force  at  each  timestep. 
By  contrast,  the  HPCMD  simulation  was  carried  out  over  25  to  64  nodes  of  an  IBM  SP2  over  the 
course  of  one  week.  The  diffusion  constant  of  0.32  ±  0.05  A^/ps  as  calculated  by  HPCMD  was 
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somewhat  lower  than  the  experimental  value  of  0.4  kVps,  but  this  deviation  may  well  be  attributed 
to  minor  inaccuracies  in  the  pair  potential  and  not  necessarily  inaccuracies  in  the  CMD 
methodology.  The  classical  system  was  a  solid  at  this  temperature.  The  preliminary  CMD 
simulations  of  liquid  hydrogen  systems  have  helped  to  establish  the  validity  of  both  the  theoretical 
and  computational  CMD  methodology.2h22 

The  second  approach  for  carrying  out  large-scale  CMD  simulations  involved  the  development 
of  effective  classical-like  pair  potentials  for  the  centroid  variables.20-22  in  this  case,  the  full  many- 
body  centroid  potential  // particles  was  approximated  as 

=  (16) 

i=\  j>i 


where  ,•,  )  was  some  effective  centroid  pair  potential  as  a  function  of  the  distance  between 

the  quantum  path  centroids  of  two  atomic  sites  i  andj;  i.e.,  r^-ij  =  \  |.  For  a  simple  pair  of 

particles,  the  quantum  centroid  pair  potentials  can  be  calculated  directly  by  performing  a  numerical 
path  integral  calculation;  i.e., 

exp  Vc,ij  fc,ij  ^■■■jDqi(r)Dqj{t)dir^^ij-rij)exp[-Sjp[qj(T),qjiT)]/n] 
where  5[q,(T),qj  ('r)]  and  5j^[qj('r),qy(T)]  are  the  action  functionals  for  the  interacting  and  free 
particle  limits,  respectively.  The  pair  separation  centroid  Ty  in  Equation  17  is  given  by 


(18) 


The  behavior  of  the  centroid  pair  potential  provided  insight  into  the  origin  of  the  dominant  quantum 
effects  in  the  solid  through  a  direct  comparison  with  the  classical  pair  potential  (see  Fig.  18  for  p- 


H2).22 


2.7  Quantum  Dynamical  Simulations  of  Solid  Hydrogen 

Historically,  there  has  been  much  interest,  both  experimentally  and  theoretically,  in  the  low 
temperature  properties  of  solid  hydrogen.  1546  xhis  interest  may  be  attributed  in  part  to  its  intrinsic 
value  as  one  of  the  prime  examples  of  a  quantum  solid,  and  to  its  practical  value  as  a  HEDM. 
Theoretically  and  computationally,  the  thermodynamic  properties  of  the  system  have  been  studied 
by  path  integral  simulations, 23-25,40  and  the  results  agree  very  well  with  experiment. 

A  CMD  calculation  was  carried  out  to  directly  determine  the  phonon  spectrum  of  pure  bulk 
solid  para-hydrogen  for  the  first  time.52  The  Silvera-Goldman  intermolecular  pair  potential,  which 
had  been  shown  to  yield  accurate  results  in  previous  simulations,  was  again  employed,‘i^2  §0  the  H2 
molecule  was  described  as  a  spherically  symmetric  particle  with  no  internal  degrees  of  freedom. 
This  approximation  is  appropriate,  since  at  the  temperature  studied  (4  K),  each  hydrogen  molecule 
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was  frozen  into  its  ground  vibrational  and  rotational  state.  The  ground  rotational  state  for  p— H2  is  J 
=  0  (7  is  a  good  quantum  number  even  in  the  condensed  phase  system),  which  is  spherically 
symmetric. 

The  PMC  and  PMD  simulations  gave  a  value  for  the  density  at  zero  pressure  that  compared 
very  favorably  with  the  experimental  value  (versus  the  classical  MD  result,  which  was  20%  in 
error).  It  should  be  noted  that  classical  simulations  at  the  experimental  density  gave  quite 
unphysical  results  (the  hydrogen  molecules  condensed  into  part  of  the  simulation  cell,  leaving  a 
vacuum  region).  Also,  the  hydrogen  molecules  were  very  delocalized  because  of  zero-point 
energy.  The  RMS  width  of  the  single  particle  distribution  function  was  18%  of  the  lattice 
spacing,25  which  was  also  what  PMC  and  PMD  simulations  predicted. 


Figure  18 

A  Comparison  of  the  Classical  Silvera-Goldman  Pair  Potential 
for  P-H2  with  the  Effective  Quantum  Centroid  Pair  Potential 
Obtained  Numerically  at  25  K 

The  shallowness  of  the  centroid  potential  in  Figure  18  effectively  reflects  the  zero-point 
quantization  of  the  hydrogen  molecules,  while  its  minimum  is  shifted  outward  from  the  combined 
effect  of  the  quantization  and  the  potential  anharmonicity. 

The  quantum  velocity  autocorrelation  functions  calculated  from  CMD  for  the  solid  at  4  K  and 
the  liquid  at  14  K  and  25  K  are  shown  in  Figure  19.  Note  the  oscillatory  behavior  of  the  function  at 
4  K,  which  is  characteristic  of  a  solid.  Classically,  the  zero  time  value  of  the  velocity 
autocorrelation  function  is  proportional  to  the  average  temperature.  This  is  not  the  case  quantum 
mechanically,  because  zero-point  energy  is  included  as  well,  which  is  temperature  independent.  It 
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should  be  noted  that  PIMC  was  chosen  as  the  method  for  calculating  the  centroid  force  for  the  solid 
hydrogen  CMD  simulations  because  it  proved  to  be  more  efficient  than  PIMD. 
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Figure  19 

A  Plot  of  the  CMD  Velocity  Correlation  Functions  for  p-ll2  for 
the  Solid  and  Liquid  at  Two  Temperatures 


The  Fourier  transform  of  the  velocity  auto-correlation  function  from  CMD  is  shown  in  Figure 
20.  This  function  can  be  related  to  the  peak  positions  in  the  experimental  phonon  spectrum  as 
determined  from  neutron  scattering  or  infrared  data.  The  CMD  calculations  agreed  remarkably  well 
with  the  experimental  results.  The  CMD  calculation  assumed  nothing  about  the  structural  form  of 
the  lattice  and  invoked  no  approximate  functional  form  for  the  phonon  motions.  It  was  purely  the 
result  of  a  many-body  quantum  dynamical  simulation.^2 


2.8  Lithium  Impurity  Trapping  in  Solid  Hydrogen 

At  the  time  of  this  report,  work  is  underway  to  study  the  actual  quantum  dynamics  of  the 
recombination  of  two  lithium  impurity  atoms  trapped  in  bulk  solid  hydrogen  using  CMD.  There  is 
a  large  thermodynamic  driving  force  for  the  separated  atoms  to  form  the  Li-Li  dimer  species,  but 
the  atoms  remain  separated  because  of  the  rigidity  of  the  intervening  hydrogen  molecules.  As  with 
other  low  temperature  hydrogen  systems,  a  classical  mechanical  treatment  of  the  problem  is 
completely  inadequate. 
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Figure  20 

The  Phonon  Density  of  States  for  Solid  p-^2  at  4  K 

The  density  of  states  in  Figure  20  was  calculated  from  the  Fourier  transform  of  the  CMD 
velocity  correlation  function  from  Figure  19. 

The  (preliminary)  centroid  potential  of  mean  force  (CPMF)  as  a  function  of  the  relative 
separation  of  two  lithium  atoms  in  solid  para-hydrogen  at  4  K  is  shown  in  Figure  21.  This 
quantity  was  computed  from  a  constant  volume  PIMD  simulation.  The  initial  configuration  was 
equilibrated  with  constant  pressure  PIMD  (zero  external  pressure)  with  the  lithium  atoms  well 
separated  in  the  lattice  in  single  substitutional  defects  (top  panel  in  Fig.  22).  In  these  defects,  the 
lithium  atoms  had  equilibrium  positions  at  the  hydrogen  lattice  positions.  The  surrounding  lattice  of 
hydrogen  was  distorted  because  of  the  larger  size  of  a  lithium  atom  (relative  to  the  hydrogen 
molecules),  but  it  still  maintained  the  basic  hexagonal  close  packed  (hep)  arrangement.  Once  the 
system  was  equilibrated,  the  volume  was  fixed  to  the  equilibrium  volume  at  zero  external 
pressure,  and  the  lithium-hthium  centroid  distance  was  holonomically  constrained.  The  mean  force 
on  the  hthium  atoms  along  the  lithium-lithium  axis  was  then  averaged  for  75  ps.  This  calculation 
was  done  for  the  entire  range  of  lithium-lithium  separations  allowed  within  the  finite  simulation 
volume  (see  lower  panel  of  Fig.  22  for  the  barrier  configuration).  The  integrated  result  of  these 
calculations  in  Figure  21  (i.e.,  the  CPMF)  gives  a  direct  estimate  of  the  quantum  barrier  for  the 
recombination  of  the  lithium  atoms.  This  barrier  was  rather  high  (=  130  K).  In  fact,  the  barrier  was 
over  30  times  the  thermal  energy  for  the  system,  which  means  the  lithium  atom  recombination  was 
a  rare  (or  activated)  event.  The  height  of  this  calculated  barrier  seems  encouraging  for  the 
possibility  of  stabilizing  lithium  impurities  in  sohd  hydrogen. 

At  the  time  of  this  report,  a  dynamical  estimate  of  the  prefactor  for  the  quantum  transition 
state  theory  (TST)  rate  constant  is  being  computed  by  CMD.  In  this  calculation,  the  lithium  atoms 
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are  equilibrated  at  a  relative  separation  such  that  they  are  in  the  well  of  the  CPMF  at  12.8  A  apart. 
The  period  of  motion  for  the  relative  lithium-lithium  distance  can  then  be  determined  from  the 
CMD  trajectories. 


Figure  21 

The  Potential  of  Mean  Force  Between  Two  Atomic  Lithium 
Impurities  in  the  Solid  P-H2  Lattice  at  4  K 


The  error  bars  in  Figure  21  are  on  the  order  of  10  K.  The  barrier  height  is  approximately  130 

K. 
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Figure  22 

Two  Snapshots  from  the  Quantum  Simulations  of  Two  Atomic 
Lithium  Impurities  in  the  Solid  p-H2  Lattice  at  4  K 


The  top  panel  in  Figure  22  is  a  configuration  from  the  stable  well  region  of  Figure  20,  while 
the  lower  panel  is  a  configuration  at  the  top  of  the  barrier. 
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3.0  WORK  IN  PROGRESS 


The  quantum  dynamical  CMD  method  is  presently  being  employed  to  computationally  probe  the 
explicit  dynamical  behavior  of  impurity  species  in  solid  para-hydrogen  through  very  largescale 
computations.  This  effort  is  focused  for  the  first  time  on  the  stability  issue  for  cryogenic  HEDM  by 
explicitly  studying  the  self-diffusion  and  recombination  steps  of  the  impurity  atoms  in  the  solid 
hydrogen  host.  The  goal  is  to  predict  the  maximum  concentration  of  impurities  that  can  be  trapped 
as  well  as  the  effects  of  temperature,  pressure,  other  impurities,  etc.  on  the  stability  of  the  material. 
In  the  broadest  context,  the  research  effort  addresses  the  critical  "scale-up"  issue  for  cryogenic 
HEDM. 


4.0  RECOMMENDATIONS  AND  CONCLUSIONS  FOR  THE  HEDM  PROJECT 

Given  the  encouraging  results  reported  in  this  document  and  the  outlook  for  future  large-scale 
quantum  dynamical  simulations  of  low  temperature  HEDM  systems,  it  is  recommended  that  the  Air 
Force  should  not  abandon  the  low  temperature  solid  hydrogen  component  of  the  HEDM  Program. 
Scientists  are  only  beginning  to  understand  the  nature  of  impurity  trapping  in  these  systems,  and  a 
wide  range  of  different  systems  and  experimental  conditions  remains  to  be  explored.  Computer 
modeling  is  poised  to  aid  in  the  search  for  viable  forms  of  low  temperature  HEDM. 
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